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Recent theoretical advances predict the existence, deep into the glass phase, of a novel phase 
transition, the so-called Gardner transition. This transition is associated with the emergence of 
a complex free energy landscape composed of many marginally stable sub-basins within a glass 
metabasin. In this study, we explore several methods to detect numerically the Gardner transition 
in a simple structural glass former, the infinite-range Mari-Kurchan model. The transition point 
is robustly located from three independent approaches; (i) the divergence of the characteristic 
relaxation time, (ii) the divergence of the caging susceptibility, and (iii) the abnormal tail in the 
probability distribution function of cage order parameters. We show that the numerical results are 
fully consistent with the theoretical expectation. The methods we propose may also be generalized 
to more realistic numerical models as well as to experimental systems. 


I. INTRODUCTION 

Upon compressions that are sufficiently rapid to avoid 
crystallization, a fluid of hard spheres (HS) first turns 
sluggish and then forms a glass [1, 2]. This glass can 
then be further compressed until the system jams [3], 
which occurs under the application of an infinite con¬ 
fining pressure [4, 5]. Glass formation is entropic, i.e., 
particles vibrate and thus cage each other in place, while 
jamming is mechanical, i.e., no motion is possible and 
particles are held steady through direct contacts with 
each other. Over the last decade, this two-transition sce¬ 
nario has been broadly validated, both numerically and 
theoretically [5-13]. Interestingly, recent advances pre¬ 
dict that - at least in the mean-field, infinite-dimensional 
{d —I oo) limit - there exists a third transition, a so-called 
Gardner transition, that is intermediate in density and 
pressure between glass formation and jamming [14-17]. 
First discovered in spin-glass models [18-22], the Gard¬ 
ner transition corresponds to a single glass metabasin 
splitting into a complex hierarchy of marginally stable 
sub-basins. The transition is thus akin to the spin-glass 
transition of the Sherrington-Kirkpatrick (SK) model, 
wherein a critical temperature separates a paramagnetic 
phase, in which a single thermodynamic state exists, from 
a marginal phase, in which a large number of distinct 
spin-glass states appear [23]. In structural glasses, how¬ 
ever, the high-temperature phase corresponds to a given 
glass metabasin that has been dynamically selected by a 
quenching protocol; it is this metabasin that then under¬ 
goes a spin-glass-like transition [24]. 

The discovery of a Gardner transition in glasses has al¬ 
ready markedly advanced our theoretical understanding 
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of jamming by providing analytical predictions for the 
critical jamming exponents [15, 16, 25-27]. It further 
suggests an explanation for the abundance of soft vibra¬ 
tional modes in glasses [27, 28], for the peculiar behavior 
of the specific heat in quantum glasses [14, 25], and var¬ 
ious other transport and thermodynamic properties in 
this regime. 

Before these fascinating problems can be tackled, how¬ 
ever, a crucial question is whether the Gardner transition 
itself, whose existence is well established in the d oo 
limit, exists in finite (low) dimensions. Renormalization 
group results indicate that the transition might disap¬ 
pear or dramatically change of nature in low d [29]. Yet 
a similar line of inquiry has been pursued for decades in 
the context of spin glasses [30], leading to the conclusion 
that, whatever the ultimate fate of the phase transition 
in the thermodynamic limit may be, the d ^ oo sce¬ 
nario provides a very good description of the system over 
the relevant experimental length and time scales [31]. At 
present, the most direct way to assess the relevance of the 
Gardner transition for the description of experimental 
glasses is through numerical simulations. It is therefore 
important to first identify the observable consequences of 
this transition in well-controlled model systems. 

This study primarily aims to develop procedures and to 
identify observables in order to reliably detect the Gard¬ 
ner transition. To that effect, we consider a simple struc¬ 
tural glass former, the infinite-range Mari-Kurchan (MK) 
model [32, 33]. The model is quite abstract and in some 
ways far from realistic models of glasses, but (i) it is a 
mean-field model by construction; (ii) it shares, in any 
finite dimension d, the same qualitative phase diagram as 
infinite-dimensional hard spheres, provided one neglects 
the effect of hopping on the glassy dynamics [34]; (iii) 
it can be studied analytically in great detail, using the 
methods of [17, 34], that we further developed for this 
work; (iv) and, most importantly, it can be easily simu- 
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lated in any finite dimensions d, including d = 3 as we 
do here. Discerning the signatures of the Gardner transi¬ 
tion in this well-controlled setting, where we are certain 
that the transition exists, shall later on enable us and 
others to study more realistic glass models for which the 
existence of the transition is not a priori guaranteed. 

The plan of this paper is as follows. In Sect. II we 
describe the MK model and its glassy behavior, and in 
Sec. Ill we detail the numerical procedures we use for the 
study. In Sec. IV we discuss several quantities that bear 
the signature of the transition, as suggested by the anal¬ 
ogy with spin glasses. For instance, at both the spin-glass 
and the Gardner transitions, the “spin-glass susceptibil¬ 
ity” diverges and the distribution of overlaps (distances) 
between different replicas becomes non-trivial. Bringing 
the various estimates of the transition together in Sec. V 
reveals that the Gardner transition can be reliably and 
reproducibly located through numerical methods in the 
MK model, and that the results are fully consistent with 
theoretical expectations. We conclude in Sec. VI with a 
description of other possible measurements to detect and 
characterize the Gardner transition, which may be more 
appropriate for numerical simulations of more realistic 
model glass formers as well as for experiments. 

Before embarking on this program, let us make a note 
of warning to the reader. Because the aim of work is to 
identify numerical methods to detect the Gardner tran¬ 
sition, we have attempted to make the numerical part as 
self-contained as possible. For what concerns the theoret¬ 
ical part, however, we have chosen to be more succinct 
and have instead relied on previous work both on spin 
glasses and structural glasses, which is here only briefly 
recalled. We expect that the reader unfamiliar with spin- 
glass theory will nonetheless be able to read and under¬ 
stand the numerical part without much difficulty. 

II. MODEL AND BASIC PHYSICAL PICTURE 

We consider a simple glass-former, the infinite-range 
Mari-Kurchan model [32, 33] - initially proposed by 
Kraichnan [35] - in which N hard spheres of diameter 
a interact through a pair potential that is a function of 
distance shifted by a quenched random vector A^. The 
total interaction energy is thus 

N 

U = ( 1 ) 

i<j 

where for particles at positions {r^} the shifted distance 
Vij is defined as = r —r^ +Aij, and u{r) is the HS po¬ 
tential, i.e., = 9{r — a) with r = jrj. The random 

shifts, which are uniformly distributed over the system 
volume V, induce a quenched disorder that suppresses 
both crystallization and nucleation between metastable 
glassy states [32-34]. The model further enables planting, 
which is a simple process for generating equilibrated liq¬ 
uid configurations at all densities [34, 36] (see Sec. Ill A). 
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FIG. 1. (Color online). Phase diagram of the HS and MK 
models in the limit d —>■ oo. Results are partially derived from 
Ref. [17], complemented by new results obtained for the Gard¬ 
ner phase [37]. The liquid EOS given by Eq. (6) (full black 
line) gives p/d = (2‘*(p/d)/2 in the limit doo. The black dot 
denotes the dynamical glass transition at 2‘^ipd/d = 4.8. Glass 
EOSs from state following (SF) for 2‘*(po/d = 5,6,6.667,7,8 
are also reported (thin colored lines, from left to right). The 
envelope for the Gardner transition for each SF (colored dot) 
is given by the dashed line. 

In all spatial dimensions, the MK model has a mean- 
field structure by construction, due to the infinite-range 
random shifts [32], and exhibits a jamming transition in 
the same universality class as standard HS [33]. The 
MK model is also fully equivalent to standard HS in the 
limit d —>■ oo, where both models can be solved exactly 
using the mean-field methods described in Refs. [5, 14- 
17, 32, 37]. In the rest of this section, we briefly de¬ 
scribe the phase diagram in the infinite-dimensional limit 
(Fig. 1), we explain how it can be applied to the MK 
model in d = 3, and we discuss some of the finite¬ 
dimensional corrections that have thus far been consid¬ 
ered [34]. 

A. Equilibrium states (liquid phase) 

The liquid phase of the MK model ergodically sam¬ 
ples equilibrium conhgurations following the Gibbs dis¬ 
tribution and has a remarkably simple structure. Its pair 
correlation function is given by 


*<■■> = ( 2 ) 

= =d(r-fT), 

where (3 = 1/T is the inverse temperature, (■ • •) de¬ 
notes thermal averaging, and denotes averaging over 
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quenched disorder, i.e., over A^. The second virial coef¬ 
ficient is 

B 2 = J j f{ri2)dvidv2 = (3) 

where the Mayer function /(r) = — 1, and Vd is 

the volume of a d-dimensional ball of unit radius. 

Because no indirect correlations exist, higher-order 
correlation functions can be factorized in a trivial way 
and the corresponding virial coefiicients are zero. For 
example, the three-body correlation function 

^ - ''>) 

= g2(r)52(r'), 

( 4 ) 

and the third virial coefficient 

B3. = JJJ .f{ri2)f{r 13)f{r23)dridr2dr3 

=-^ JJJ .f{ri 2 )f{ri 3 )x ( 5 ) 

X /(ri 3 ~ ri2 -l- Ai 2 -I- A23 — Ai3)(iri(ir2(ir3 
= 0 . 

Note that if |ri 2 |, |ri 3 | < ct, then /(ri 3 - ri 2 -I- A 12 -I- 
A 23 — A 13 ) = 0 in the thermodynamic limit because ran¬ 
dom shifts are uncorrelated and typically of the system 
size, and thus |ri 3 - ri 2 -I- A 12 -I- A 23 - A 13 I > a. It is 
straightforward to generalize this argument to show that 
all higher-order virial coefficients are also zero [32]. 

Because only the second virial coefficient is non-zero, 
the reduced pressure p equation of state (EOS) for the 
liquid is 

p = pP/p=l + B2P=l + 2‘^-^p, ( 6 ) 

where the combination of inverse temperature jd, pressure 
P, and number density p = N/V gives a unitless quantity 
p whose only dependence is on the liquid volume fraction 
(p = pVd{al2)‘^. 

B. Dynamical glass transition 

Although the structure and thermodynamics of the liq¬ 
uid are trivial, its dynamics is not. We will focus here 
on the equilibrium dynamics, i.e. starting from an equi¬ 
librium initial condition [38]. In infinite dimension, a 
dynamical glass transition (pd separates two distinct dy¬ 
namical regimes. For p < pd^ the dynamics is diffusive at 
long times, as expected of any liquid. Upon approaching 
Pd, however, the dynamics grows increasingly sluggish, 
and above pd, each particle is fully confined within a cage 
formed by its neighbors. The typical size of that cage is 
the cage order parameter Ai (the meaning of the suffix 


will become clear below), which, in that regime, can be 
extracted from the long-time limit of the mean-squared 
displacement (MSD) 

1 ^ _ 

(7) 

Ai = lim A(t) . 

i—>-oo, 

From the d —>■ 00 solution, we know that the equilib¬ 
rium distribution of the order parameter, Peq(A), has two 
peaks for p > pd (see Fig. 2) [39, 40]. The first character¬ 
izes the distance between two glass configurations within 
a same metabasin. It is centered around Ai, which is the 
typical size of this basin. The second characterizes the 
inter-basin distance. It is centered around Aq = 00 , be¬ 
cause states that belong to different metabasins are com¬ 
pletely uncorrelated. In technical terms, this situation is 
described by a 1-step replica symmetry breaking (IRSB) 
scheme [23]. Note, however, that the peak at Ai has an 
exponentially small weight in N, because there exists an 
exponentially large number of distinct glass states [40]. 
Hence, in the thermodynamic limit, Peq(A) = 6{A — Aq) 
everywhere in the liquid phase, i.e. even for p > pd- 

C. Glass state following and the Gardner transition 

In d —>■ 00 , each equilibrium configuration at density 
‘Po > Pd is forever trapped into one of the exponentially 
many glass metabasins. The pressure of an equilibrium 
configuration at po is given by the liquid EOS, Eq. (6), 
but if one compresses (or decompresses) such a configu¬ 
ration up (or down) to a density p^ the system remains 
within the metabasin that was initially selected. The 
system thus falls out of equilibrium in the sense that it 
cannot visit the ensemble of all possible distinct glass 
metabasins and thus follows an EOS different from that 
of the liquid. 

The infinite lifetime of these d —)■ 00 glass states im¬ 
plies that one can adiabatically follow the EOS of a single 
glass state through a restricted equilibrium approach, a 
construction also known as state following (SF) [17, 41- 
43]. SF consists in using a first configuration that is 
equilibrated at initial density i^oi and a second config¬ 
uration that is in a restricted equilibrium at a different 
density p [17, 41-43]. The restriction is that that sec¬ 
ond configuration must be part of the same metabasin 
as the first one. The state of the second configuration 
thus describes the evolution of the glass metabasin with 
p. Dynamically, this process corresponds to preparing 
a system at equilibrium at initial density po, then com¬ 
pressing it at density p, and assuming that it is able to 
equilibrate inside the glass metabasin, but without escap¬ 
ing it. In simpler words, structural relaxations are frozen 
and particles only rattle inside their cages. Upon com¬ 
pression, these cages shrink until the jamming density, 
Pj(Pq), is reached. Particles are then mechanically in 
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similar to the ones reported in Ref. [16], but the full de¬ 
tails of this work are given elsewhere [37]. 

Following a glass state in restricted equilibrium gives 
for (f < (pG{(po) a distribution Psf(A) that has a single 
peak at Ai. Note that because the system is confined 
to a single glass basin, the peak at Aq is absent. At 
‘Pci^Po) the glass basin fractures in a SF-fullRSB struc¬ 
ture, and correspondingly the single peak in the distri¬ 
bution Psf(A) fractures into two peaks centered around 
Aea and Ai, connected by a wide continuous band (see 
Fig. 2 and [23]). Here, Aea is the typical size of the inner¬ 
most sub-basins at the lowest hierarchical level, and Ai 
is the typical distance between the outermost sub-basins, 
i.e., the size of the metabasin. For the same reason as 
above, Psf(A) does not show the Aq peak. 

Because the cage order parameter changes continu¬ 
ously at ipG, the Gardner transition is a continuous 
critical transition [18, 19]. The Gardner phase is also 
marginally stable, in the sense that a zero mode is always 
present in the stability matrix of the free energy [23]. Be¬ 
cause jamming is located within the Gardner phase, its 
marginal stability and critical scaling behaviors can con¬ 
sequently be obtained from a fullRSB thermodynamic 
calculation [15, 16]. 


FIG. 2. (Color online), (a) Organization of glass states (blue 
dots) for yjQ > Pd', the typical intra-basin MSD is Ai, while 
the inter-basin MSD is Aq = oo. (b) In the SF-fullRSB phase, 
pa < p < Pj , glass metabasins subdivide into a hierarchical 
structure of sub-basins with typical innermost MSD Aea and 
outermost (meta)basin MSD Ai. Schematics of the equilib¬ 
rium Peq(A) and restricted equilibrium Psf(A) (blue area) 
distributions are given for (c) the SF-RS (or IRSB) and (d) 
the fullRSB phases. 

contact, which makes the system mechanically rigid. In 
the following, we call the restricted equilibrium EOS of 
a given glass metabasin simply its glass EOS. 

Examples of different d ^ oo glass EOS obtained by 
SF are given in Fig. 1 [17, 37]. In d —>• oo, a compressed 
state under restricted equilibrium undergoes a Gardner 
transition at pg(pq), at which point the glass metabasin 
subdivides into a hierarchy of sub-basins. Technically, 
before the Gardner transition, i.e., for p < pg{pq), 
the glass metabasin is obtained by a replica symmet¬ 
ric SF (SF-RS) computation [44], while in the Gardner 
phase, i.e., for p > pg{pq), a full replica symmetry 
breaking SF (SF-fullRSB) computation is needed [15- 
19, 23]. The Gardner transition is therefore akin to 
the spin glass transition one of the SK model, but re¬ 
stricted to a single glass metabasin. The high tempera¬ 
ture phase of the SK model thus corresponds to a sim¬ 
ple glass metabasin for p < pg(po) while the low tem¬ 
perature phase corresponds to a fractured and marginal 
metabasin for p > pg(po)- Note that we here com¬ 
plement the SF-RS computation of Ref. [17] with a SF- 
fullRSB computation, in order to obtain the complete 
EOS reported in Fig. 1. The SF-fullRSB equations are 


D. Timescales 

The timescales that characterize the dynamics in the 
different phases of the MK model are predicted based on 
the general correspondence between statics and dynam¬ 
ics in spin glasses [45, 46]. Here again, we only consider 
the system behavior in the limit d ^ oo. Beyond the mi¬ 
croscopic timescale tq, over which dynamics is essentially 
ballistic, one gets the following picture. 

• In the liquid phase below p^., dynamics is character¬ 
ized by two timescales: the /3-relaxation timescale 
Tp, over which particles explore their transient cage, 
and a longer timescale Tq,, over which dynamics is 
diffusive. Both timescales are finite for p < p^., and 
diverge at pd according to the scaling predicted by 
mode-coupling theory (MGT) [47]. We do not dis¬ 
cuss this regime further because it is not directly 
related to the Gardner transition. 

• In the liquid phase above pd and in the SF-RS 
(simple glass) phase, the same two timescales exist: 
Tq ~ exp(A) is the timescale for jumping from one 
glass metabasin to another, which being infinite in 
the thermodynamic limit properly defines the glass 
metabasins, and Tp is the timescale for equilibrat¬ 
ing within a glass metabasin, i.e., for particles to 
explore their cage, which is finite for p > pd but 
diverges upon approaching pd from densities above 
it, again following a MCT-like scaling form [47]. 

• Upon approaching the Gardner transition, rg again 
diverges, scaling as Tp ~ {pG — pV■ It does so 
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for the exact same reason as at near the spin-glass 
transition of the SK model and corresponds to the 
critical slowing down close to a second-order phase 
transition. 

• In the Gardner phase, dynamics is described by 
three timescales: Tp is the timescale for equilibrat¬ 
ing within a single glass sub-basin (A < Aea)) 
Tmeta is the timescale over which system explores 
the structure of the sub-basins within a given glass 
metabasin (Aea < A < Ai), and ~ exp(A) 
remains the timescale for jumping from one glass 
metabasin to another (A ~ Aq). To be more clear, 
none of these processes correspond to a simple ex¬ 
ponential with a single timescale. Everywhere in 
the Gardner phase Tp = oo, because of the phase’s 
marginality. The relaxation inside a single sub¬ 
basin is thus expected to scale as a power-law in 
time. The exploration of sub-basins is character¬ 
ized by a complex distribution of free energy barri¬ 
ers and relaxation times, and thus Tmeta ~ exp(A“) 
with Of < 1 (a = 1/3 is expected in the SK 
model [23]). Note that barriers between sub-basins 
are much lower than those between metabasins, 
hence Tmeta < Ta. 

Let us now situate the scheme for our work in this 
complex d —>■ oo dynamic phase diagram. We are 
conveniently not concerned with the scalings of the 
MGT regime. Thanks to the planting procedure (see 
Sec. Ill A), we are indeed able to easily generate equi¬ 
librium configurations of the MK model for arbitrary (/?o, 
including lpq > These configurations are well equi¬ 
librated in a glass metabasins and therefore have an in¬ 
finite Tq, (in the thermodynamic limit) and a finite t^. 
Because = oo in this density regime, we can simply 
forget about its existence and consider that we are for¬ 
ever restricted into a given glass metabasin. From this 
point of view, we are in a situation similar to that of a 
spin glass where one is able to start at equilibrium in the 
paramagnetic, high temperature phase. 

If we compress the system slowly enough to a hnal 
density ip G {(po, (po), then equilibration within the glass 
metabasin is possible, because remains finite. In this 
situation, the system behavior is described by the SF- 
RS computation. To study the Gardner transition, how¬ 
ever, we need to compress the system up to ipa, at which 
T /3 diverges. Hence, upon approaching tpo, Tp eventu¬ 
ally becomes larger than the simulation timescale, and 
equilibration (even in the restricted SF sense) becomes 
impossible. The system thus falls out of (restricted) equi¬ 
librium and is not described anymore by the SF compu¬ 
tation. 

For p> > {pQ the situation is even worse both as Tp and 
Tmeta are infinite. The SF-fullRSB computation, which 
gives the restricted equilibrium properties in this regime, 
is only an approximation even at long (but not divergent 
with exp(A)) times. The situation is akin to that in 
the SK model where the fullRSB computation is only an 


approximation of the states reached dynamically at long 
times in the spin-glass phase. Because in this regime 
the planting technique does not work (it only works on 
the liquid line in Fig. 1), we cannot use it to study the 
restricted equilibrium in the Gardner phase. 

In the following we will therefore present two kinds of 
data: 

• for (^0 A far enough from ipQ such that Tp 

is smaller than the simulation timescale, we obtain 
restricted equilibrium data 

• for ~ ipQ and p) > ipc, the system is out of 
(restricted) equilibrium, shows aging effects, and, 
at long times, we obtain states that are qualita¬ 
tively similar to the SF-fullRSB ones, but not ex¬ 
actly equal to them. 

In Sec. Ill we give a more precise definition of the nu¬ 
merical protocol we use to study these different regimes. 
Note that, according to the above discussion, from now 
on we will refer to the restricted equilibrium simply as 
“equilibrium”, given that we always work inside a glass 
metabasin. 


E. Finite-dimensional MK model 

To conclude this section, we discuss the additional ef¬ 
fects that appear when one considers the MK model in 
finite dimensions, and that affect the above discussion. 

From a theoretical point of view, a hnite d affects quan¬ 
titatively the phase diagram, but to a lowest degree of 
approximation one can simply take the phase diagram of 
Fig. 1 and fix the value of d (e.g. d = 3) to obtain re¬ 
sult for the EOS. However, there are several systematic 
corrections that impact the accuracy of this result. 

• The infinite-dimensional results are obtained 
within a Gaussian structure of the cage, which pro¬ 
vides the exact results in the limit d —)■ oo [48]. 
The Gaussian equations are slightly different in fi¬ 
nite d [5]. The corrections have the form of a se¬ 
ries in 1/d and are quite small. For example, in 
d = 3 one has 2^(/?d/3 « 4.74 instead of the infinite¬ 
dimensional result 2'^(/?d/d « 4.8. These corrections 
are negligible and could be easily taken into account 
if needed. 

• A more important problem is the inexactitude of 
the Gaussian assumption for the MK model in finite 
dimensions. One should instead optimize the free 
energy over a generic cage function, but this com¬ 
putation is technically quite difficult. In Ref. [34], 
however, two different cage functions were studied 
and found to give similar results. The corrections 
coming from the non-Gaussianity of the cage are 
thus also rather small, although they are more dif- 
hcult to estimate. 
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• In finite d, cages are heterogeneous [34]. Care¬ 
fully taking this effect into account would require a 
cavity calculation similar to the one performed in 
Ref. [34], which is beyond the scope of this work, 
but is also likely a small correction for tpo > [34]. 

For all these reasons, the results of Fig. 1 are only an 
approximation to the finite-dimensional phase diagram. 
In the following, we will nonetheless take the results of 
the simplest theoretical approximation, namely take the 
d —)■ oo results of Fig. I and use them in d = 3, after 
properly rescaling the axes. We expect (and check) that 
the corrections are quite small. Sometimes, however, in 
order to take effective account of these corrections, we 
will further rescale the results, so as to obtain a better 
quantitative agreement between theory and simulation. 

The most important genuinely finite-dimensional ef¬ 
fect is hopping. As discussed in Ref. [34], slightly above 
iy9d, particles are not perfectly confined in their cages, 
as one would expect based on the d —>■ oo picture. In¬ 
stead, each particle is allowed to explore a network of 
cages, connected by narrow pathways. Hopping consists 
in particles jumping between distinct cages. The corre¬ 
sponding time scale is finite at (/ 9 d and thus the d —>■ oo 
divergence of the relaxation time is washed away in fi¬ 
nite d. Hopping effects also change deeply the dynamics 
of the system with respect to the d —)■ oo prediction, 
and the configurations at (^o ^ Vd are not able to con¬ 
strain the dynamics for infinitely-long times. Because the 
timescale for hopping increases quickly upon increasing 
both density and dimension, considering values of (^o that 
are slightly above i^d suffices. In practice, based on the 
analysis of Ref. [34], in d = 3 and for (/?o > 2.5 hopping is 
strongly suppressed, hence in that regime the mean-field, 
d —)■ oo scenario should apply reasonably well. 


III. NUMERICAL APPROACH 

In this section, we provide the numerical details used 
in the simulations of the glass states of the MK model. 


A. Planting 

An important algorithmic advantage of the MK model 
is that planting can be used to generate equilibrium liq¬ 
uid configurations at any lpq [34, 36]. This procedure 
sidesteps the tedious and time-consuming work of first 
preparing dense equilibrium configurations, as would be 
needed for typical glass formers, such as HS. The basic 
idea is to switch the order in which initial particle posi¬ 
tions {ri} and random shifts {A^} are determined. One 
first chooses the particle positions {ri} independently 
and uniformly in the volume U, and then for each particle 
pair one chooses a random shift A^ uniformly under the 
sole constraint that the two particles should not overlap, 
which is quite straightforward to satisfy. As long as the 


quenched and the annealed averages of the free energy are 
the same (see Ref. [36] for a more detailed discussion), 
a planted state is a true equilibrium state and automat¬ 
ically satisfies the liquid EOS, Eq. ( 6 ). This condition 
is met along the replica symmetric phase for < ipY.-, 
where i^k is the Kauzmann point at which the config¬ 
urational entropy vanishes [5, 32]. Because in the MK 
model (/?K = oo [32], planting a liquid configuration is 
thus possible at any density, which dramatically reduces 
the computational cost of the initial equilibration. 

In our notation, a given {r^} and {A^} defines a sam¬ 
ple. A sample thus identifies a given system (defined by 
{Ay}) and, for this system, one of its glass metabasins 
(selected by {ri}). 


B. Molecular dynamics (MD) simulations 

We use event-driven molecular dynamics (MD) to sim¬ 
ulate MK particles in d = 3 [7, 34]. Periodic boundary 
conditions with the minimum image convention are im¬ 
plemented on the shifted distances [r^ — Vj -|- Ay |. Time t 
is expressed in units of \/where the particle mass 
m and diameter a as well as the inverse temperature /? are 
set to unity. Systems consist of A = 800 particles unless 
otherwise specified. This system size is large enough to 
contain a first full shell of neighbors around each particle, 
and to keep the periodic boundary effects on caging to 
a minimum [34]. Einite-size effects are studied for some 
of the observables for the initial liquid density po = 2.5 
(see Sec. IV B). 

To simulate SE, for a given sample, we start from the 
planted equilibrium configuration at a packing fraction 
ipo, and grow the spheres following the Lubachevsky- 
Stillinger algorithm [ 6 , 7] at constant growth rate 7 = 
0.001, unless otherwise specified, up to a desired p. Once 
compression is stopped at the target density, the origin 
of time is set. We will thus typically (although not al¬ 
ways) define the waiting time tyj, as the time that has 
elapsed since the end of the compression. Erom that mo¬ 
ment on we start measuring observables, keeping density 
and temperature (and thus energy) constant. Note that 
because 7 is finite and rather small, part of the equilibra¬ 
tion happens already during compression, so provided we 
are not too close to po the system is stationary at all tw 
(we come back on this point later). This procedure is re¬ 
peated over Ns samples in order to average over thermal 
and quenched disorders. Errors are computed using the 
jack-knife method [49]. Depending on the statistical con¬ 
vergence of the different observables, Ns is varied from 
500 to 75,000, as specified in the discussion of the various 
measurements. 
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C. Observables 


The pressure evolution along SF is reasonably well de¬ 
scribed by a free-volume EOS 


1 

P 


= C 



( 8 ) 


where C is a fitting parameter. Fitting Eq. ( 8 ) to the 
compression results provides an estimate of (/ 7 j (see Ta¬ 
ble II). If a sufficiently small 7 is chosen, no aging is 
observed in the pressure, and using slower compression 
rates gives only negligible corrections to the glass EOS 
(see Fig. 3). Interestingly, upon decompression, the state 
follows the same EOS up to a threshold density at which 
it melts into a liquid phase. This phenomenon has been 
recently predicted by the theory [17, 50], and observed 
numerically in simulated ultrastable glasses [51, 52]. 

To obtain more structural information about the free 
energy landscape, we also simulate a cloning procedure. 
The approach consists of taking two exact copies (clones) 
A and B of the same planted configuration at tpo, and as¬ 
signing them different initial velocities, randomly drawn 
from the Maxwell-Boltzmann distribution. These two 
copies are then independently compressed up to (p, be¬ 
fore measuring the mean-squared distance between them 



FIG. 3. (Color online), (a) Compressions (7 > 0) and decom¬ 
pressions (7 < 0) of an initial equilibrium state at po = 2.5 
and 4.0. The results are averaged over Ng = 100 samples. 
The theoretical curves uses data from Ref. [17] and from this 
work [37]. (b) Compression data are indistinguishable for 

7 < 0 . 01 . 


IV. DETECTING THE GARDNER 
TRANSITION 


1 ^_ 

= ( 9 ) 

where rf{t) and rf (t) are the positions of particle i at 
time t in clones A and B, respectively. Although the 
two clones start from the same initial configuration, their 
compression histories are different once p is reached. 

The detailed behavior of the clones will be discussed in 
Sec. IV A, but let us explain here briefly why the cloning 
procedure is useful to detect the Gardner transition. In 
the SF-RS phase, the two clones are uncorrelated in the 
glass basin and Aab {t) converges quickly (on a time scale 
t ~ if the two clones are not sufficiently well equili¬ 
brated along the compression) to the equilibrium value 
Aab = Ai. Hence Aab = linit-j-oo A(t) in the SF-RS 
phase. By contrast, if the end point of the compression 
falls within the SF-fullRSB (or Gardner) phase, clones 
most likely fall into different sub-basins. Their mean 
square distance can then be described by a non-trivial 
time-dependent probability distribution, PABit, A), that 
depends on the way sub-basins are sampled. Calculating 
these weights is difficult, because the two clones are gen¬ 
erally out of equilibrium. Because the probability that 
the two clones fall in the same state is very small, how¬ 
ever, and therefore AAB{t) > Aea for all t. Hence, in 
the Gardner phase Aab)!) at all t is strictly larger than 
the long-time limit of the MSD. The long time limit of 
Aab ~ A, being zero in the SF-RS phase and non-zero in 
the SF-fullRSB phase, thus provides an order parameter 
for the Gardner transition. 


In this section, we describe different means of detect¬ 
ing the Gardner transition through numerical simulations 
designed to follow the evolution of glassy states. Let us 
stress once again that the aim of this work is to under¬ 
stand, in a controlled setting, how well numerical sim¬ 
ulations and experiments can detect the Gardner tran¬ 
sition. We can therefore test different strategies to un¬ 
derstand their advantages and limitations, and use the 
analytical results to assess the quality of the numerical 
results. Based on the discussion in Sec. H, we follow 
two complementary approaches. The hrst is based on 
dynamics. From the long-time dependence of the mean 
square displacement, we determine rg, whose divergence 
in the glass metabasin signals the Gardner transition. 
The second is based on the properties of the distribu¬ 
tion function Pab(A), which becomes non-trivial in the 
Gardner phase. We investigate different moments of this 
distribution to detect signatures of the transition. 


A. Dynamics 

As discussed in Sec. HD, the Gardner transition is a 
second order phase transition associated with a diverg¬ 
ing characteristic relaxation time rg that controls the 
dynamics of the MSD. To detect this relaxation time 
we make use of the following observables already briefly 
discussed in Sec. HI. Recall that the origin of time 
t = = 0 is set at the end of the compression when 

a given final density p is reached. 









1. Definition of the relevant observables 


We define the MSD between two configurations at dif¬ 
ferent times: 

1 ^ 

+ , (10) 

' i=l 

and the MSD between two different clones 
1 ^ 

^Asit) = J^^\rf(t). ( 11 ) 

' Z =1 

From these two instantaneous quantities we can define 
different observables. The statistical average (over com¬ 
pressions and samples) gives 


A(t,t„) = , AABit) = (AAB(t)) , (12) 

and we define 


6A{t,ti,) = AAB{.t + t^) - A{t,tni) ■ (13) 


We also define a time-dependent caging susceptibility as 
the normalized variance of the MSD 


(A2(t,tu,)\ - (A(t,tu,) 

x(t,U) = ^ 2 - 

(A(t,t„)) 

and its counterpart of cloned configurations 


{A\^{t)) - (AABit) 
XABit) = ^ \ , - 

(AABit)'^ 


(14) 


(15) 


In the SF-RS phase in restricted equilibrium, i.e., at 
large enough tw and t > 0, A(t,tu,) = A{t) gives back 
Eq. (7), while A^B(t) = Ai does not depend on time. In 
this case 


Ai = AABi'^t) = lim A(t) (16) 

t^OO 

gives the average cage radius of the glass basin. Similarly, 
= x{t) and XAB{t) = X, where 

X = XABi'^t) = lim x(l) (17) 

is the average susceptibility of the glass basin. Finally, 
note that based on Eq. (16), we have 

lim lim SA{t,t^)=0. (18) 

t—^OO t-uj—^OO 

By contrast, in the SF-fullRSB phase, equilibrium is not 
reached even for very large (recall that we do not con¬ 
sider here times that are comparable to exp(A)). There¬ 
fore, dA{t,tjju) remains non-zero even for large tw and t, 
which makes the long time limit of dA{t,t-w) a dynamic 
order parameter for the Gardner transition. 



10-^ 10-1 10" lOi 10^ 


t 

FIG. 4. (Color online). SF from ipo = 2.50 gives (a) AAs(t), 
(b) A(t,tu, = 0), and (c) SA{t,tuj = 0) at different ip. Solid 
lines are fits to Eq. (21) (bottom panel). 


2. Qualitative change in caging and susceptibility 


Figure 4 shows the time dependence of AAB{t), 
A{t,tjjj = 0) and SA{t,tjjj = 0) for (po = 2.5, averaged 
over Ns = 15,000 to Ng = 500 from the lowest and high¬ 
est (p, respectively. Both A and Aab are observed to be¬ 
have slightly differently below and above pc {pc ~ 3.00 
for pq = 2.50, a more precise estimate is obtained be¬ 
low). For p < Pq, A(t,tjj, = 0) first (and up to a mi¬ 
croscopic time To ^ 10”^) grows quickly because of the 
ballistic motion of particles [34] and then more slowly in 
the relaxation regime, before eventually reaching the 
plateau A = Ai that defines the cage size. This plateau 
coincides with the (almost time-independent) results for 
Aab (t), as is qualitatively expected for a system in a SF- 
RS phase. For p > pQ, the situation is a bit more convo¬ 
luted. As discussed above, beyond the Gardner transition 
each of the original metabasins is expected to subdivide 
into a hierarchical distribution of glassy states. Equili¬ 
bration within the glass metabasin is now, however, im¬ 
possible, and we thus observe that Aab(0 depends on 
time for all observable times while A(t,tuj = 0) remains 



























9 




10-1 10 “ 10 “ 10 “ 10 “ 


t 

FIG. 5. (Color online), (a) Long-time results for 
The data show that, for ip < ipa, AAs(t) saturates at its 
asymptotic value (dotted line) in the long-time limit (see, e.g., 
<p = 2.88,2.90). (b) A(t, tiu) and (c) hA(t, t™) are plotted as 
a function of t for various and t™. 

always strictly smaller than AAsit) and never reaches a 
plateau. Correspondingly SA{t,tw = 0) does not decay 
to zero. 

In Fig. 5, we select a few densities and consider the evo¬ 
lution of AAsit) over much larger times than in Fig. 4 as 
well as the dependence of A{t, t^) on both t and Note 
that these data are averaged over many fewer samples. 
Ns « 300. We consider the detailed behavior of these 
two quantities: 

• In general, AAsit) (Fig. 5a) should correspond 
to the average distance between conhgurations re¬ 
stricted to within a given metabasin, but for t < 
Tp the basins are sampled with non-equilibrium 
weights, and hence AAsit) slowly drifts. For ip < 
ipQ, we indeed observe that once t ^ Tp, AAsit) 
reaches a stationary value. Upon approaching (pc 


and for (p > (pc, however, Tp is becomes so large 
that the drift of AAsit) persists at all simulated 
times. Note that the drift can be positive or nega¬ 
tive, depending on density. 

• The behavior of A{t,tyj) (Fig. 5b) is naturally de¬ 
scribed as a function of t for hxed tw Below (pc 
(e.g., for (p = 2.85), A{t,tyj) is independent of 
tw and behaves as in Fig. 4. Beyond (pc (e.g., 
for (p = 3.10), however, the system is initially 
trapped into a sub-basin, and thus A{t,t^) grows 
until reaching a plateau corresponding to the size 
of the sub-basin, Aea 7 for t ^ Tineta(tiu)- For 
t Tineta(tuj), the System can explore the struc¬ 
ture of sub-basins, hence A{t,tyj) keeps increas¬ 
ing, and no clear first plateau can be detected. 
A second plateau should be reached in the limit 
t >> Tineta(tiu)) when the metabasin is fully ex¬ 
plored, but this regime is here beyond computa¬ 
tional reach. The timescale r^etaitw) indeed in¬ 
creases with tw [45, 46], as is clearly visible in 
Fig. 5, where the drift shifts towards larger times 
upon increasing tw Note that in the limit N —> oo, 
Tmetaitw) should diverge for tu, —>■ oo, but a finite 
N acts as a cutoff and Tjneta{tw) should instead sat¬ 
urate to a value ^ exp(A^/^) for large 

• In the Gardner phase A{t,tw) < AAsit) for all ac¬ 
cessible times Therefore, although (5A(t,tu,) 

goes to zero at large times for (p < (pc, for (p > (pc 
it does not fully decay in the accessible t regime (for 
all tyj) and instead converges to a plateau (Fig. 5c). 
Note also that the evolution of 6A{t,t,„) with in¬ 
creasing density is qualitatively identical to that 
of the overlap with decreasing temperature in spin 
glasses, as reported in Ref. [53, Fig. 8 ] (compare 
also to Fig. 4c). These results further support the 
strong analogy between the Gardner and the spin- 
glass transitions. 

The dynamical behavior of the caging susceptibility 
is qualitatively similar that of the MSD (Fig. 6). As 
for the cage order parameters, in the SF-RS phase (for 
(p < (pc) the two susceptibilities become identical in the 
long time limit, i.e., x(t — = XAsit oo). 
By contrast, in the SF-fullRSB phase (for (p > (pc), 
on a timescale t < Tyneta{tw), we generally observe that 
x{t,tw) < XAB{t)- Note that the magnitude of the sus¬ 
ceptibility increases by more than a decade in the density 
range considered here, which is a clear signature of the 
Gardner transition. A detailed analysis of this increase 
is discussed in Sec. IV B 3. 


3. Computation of timescales 

Based on the qualitative picture of the dynamics 
presented above, we now attempt a more quantitative 
analysis based on the classic work of Ogielski on spin 
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FIG. 6 . (Color online), (a) Time evolution of the caging sus¬ 
ceptibility = 0) (scatters) and XAsit) (lines), (b) For 

three selected densities, we report the evolution of XAB{t) 
for much longer times (dotted lines) and the dependence of 
x(t,tm) on both t and t-w (see legend). 


glasses [53] and its recent extension to the study of the 
dynamical transition in the d = 3 Edwards-Anderson 
model under an external field [54]. The idea consists 
in obtaining a relaxation timescale T(tu,) from the decay 
of 6A{t,tiu) (Fig. 4c). Note, however, that this scheme is 
only well defined in the SF-RS phase, where equilibrium 
can be reached and Eq. (18) holds. In the SF-fullRSB 
phase, SA{t,tw) does not decay to zero and evolves con¬ 
tinuously over a broad range of time scales. Extracting a 
single timescale is then not so straightforward. We will 
focus in the following on the data for tw = 0 , for which 
we have more statistics. For ip < ipQ, these data are also 
representative of all tw 

In order to facilitate the numerical analysis, we make 
use of analytical results. According to the general the¬ 
ory of critical glassy dynamics developed in Refs. [55-57], 
upon approaching the Gardner point in (restricted) equi¬ 
librium from the SF-RS phase, hence when SAit^t^) = 
SA{t) does not depend on one has: 

SA{t) ~ 5pT{tjT/3) , Tp ~ Sp~'^ , (19) 

where dp = \p — and the function T{x) is such that 
J^(a: <C 1 ) ^ x~°‘ while J^(a: » 1 ) decays exponentially. 
Here the exponents a and 7 = 1/a are related to the 


2‘^Pold 

2 VG/d 

A 

7 = 1 /a 

4.8 

4.8 

0.7027 

3.069 

4.9 

5.34 

0.5607 

2.651 

5 

5.64 

0.5091 

2.547 

5.25 

6.18 

0.4378 

2.427 

5.5 

6.61 

0.3938 

2.363 

6 

7.33 

0.3398 

2.295 

6.667 

8.16 

0.2957 

2.245 

7 

8.54 

0.2801 

2.228 

8 

9.63 

0.2469 

2.194 

10.667 

12.36 

0.2042 

2.154 


TABLE I. Analytical results (in the limit d -A 00 ) for po and 
the exponents A, 7 and a for several po . 


so-called MCT exponent parameter A by the relation 


r(i-a)^ 

r(l-2a) ■ 


( 20 ) 


The parameter A can be computed analytically within the 
replica method by analizing the cubic terms of the replica 
action [14, 56, 57]. The results of this computation are 
reported in Table I (details are given in Refs. [37]). 

In order to estimate Tp, we fit the results for SA{t, = 
0 ) using an empirical form F{x) oc x~°‘e~^ that has been 
used for spin glasses [53, 54], which gives 


6A{t,tiu 


0 ) = 


exp[-(f/T^)^ 


( 21 ) 


where the parameters a, b and Tp depend on p, po, and 
offers a first estimate of Tp. Note that we fit the expo¬ 
nent a instead of using the analytical result, because the 
critical regime over which the exponent coincides with 
a is narrow and away from pc the effective exponent is 
quite different (see Ref. [53, Fig. 12]). With this choice 
all the fits are very good, as reflected by the Pearson Xp 
per degree of freedom (d.o.f.) [58] being much less than 
1 (examples of the quality of the fit are given in Fig. 4c). 
As p approaches pc for a given <po, Tp is expected to 
diverge as 


Tp ~ \p-Pc\ (22) 

In order to obtain a more constrained value of Pc^ we 
fix the exponent 7 to its analytic value given in Table I 
and only fit p'q and the prefactor. This time the fits are 
not excellent, with values of Xp/d.o.f. of the order of 1 
or larger (see Fig. 7). The results for p'q are reported in 
Table II. 

An alternative estimate of Tp can be obtained from 
the logarithmic scaling of 5A{t,tu, = 0 ) at long times 
(see Fig. 8 ). For p < pc ^ the fitting form 


5A{t,tw 


0) =fc 


log(0 

log(T-^') 


(23) 
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FIG. 7. (Color online), (a) Growth of r'p with — ip for ipQ 
from Table II. (b) The two estimates of rp, r'p and Tp, as a 
function of Pq — for ipo = 2.5. (c) Evolution of a and b with 
ip for po = 2.5. 


0.0025 


0.002 


’ 0.0015 


< 0.001 


0.0005 






* * 




if = 2.80 * 

ip = 2.85 o 
if = 2.90 • 
if = 2.95 
if = 2.97 ♦ 


i I 




* -4 








* ^ 


10 ' 


® ® m 

i I s ® S 

—— ' » se - 

10 " 


^ i 4 


FIG. 8. (Color online). Logarithmic decay of SA{t,tjn = 0) 
with time for po = 2.5. The results are fitted to Eq. (23) 
(dashed lines). 


observables unless otherwise specified. 


with a density-dependent constant k gives Tp [54]. Com¬ 
paring Tp and Tp suggests that the divergence of the two 
timescales is compatible with a same pc and a very sim¬ 
ilar power-law exponent (see Fig. 7b). The insensibility 
of the estimator of Tp to its precise definition adds sup¬ 
port to our claim that the observed divergence is due to 
a true thermodynamic transition. Also, we repeated the 
analysis for > 0 with very similar results. 

To conclude the discussion, note that the results for 
AAB{t) and A(t,tuj) at low densities, i.e., Pd < Po ^ 2.2, 
may be affected by hopping (Sec. HE and [34]). In these 
systems the timescale for leaving a metabasin (albeit only 
through local hopping processes) is comparable to Tp 
even near pQ. The estimate of po in this regime is there¬ 
fore subject to a larger error, which explains the bigger 
difference between p'^ and other pQ estimates (Table II). 
For the limit case pd = 1.8, we do not even attempt to fit 
the data because no clear power-law regime can be distin¬ 
guished. By contrast, for p^ > 2.5, hopping is negligible 
on the timescales achieved numerically. 


1. Average mean square displacement 

We consider the averages of A and Aab and com¬ 
pare them with the theoretical SF results (see Fig. 9 for 
Pd = 2.5). In order to take into account the corrections 
discussed in Sec. HE, both the numerical and the theo¬ 
retical datasets have been rescaled. For the vertical axis, 
we rescale the results to Ai((^o)j which is known both 
in theory and in simulation, such that both datasets are 
equal to 1 for p = pQ. For the horizontal axis, we rescale 
p to pG with the theoretical pQ for the SF-RS curve 
(Table I) and a fit factor Pq for the numerical results 
(Table II). Figure 9 clearly shows the difference between 
the two expected regimes (Sec. IIC): for p < pc, we 
obtain A ^ A^^ ~ Ai, while for p > pQ, we obtain 
A ^ Aea and Aab Ai, and hence Aab > A. Al¬ 
though at short times we observe AAB{t) « Ai, on much 
longer timescales we expect AAB{t) to evolve slowly to¬ 
wards its equilibrium value AAB^t) = (A)sf, which is 
however only approached for t that diverge with N. 


B. Static functions 

As we have shown in Sec. IV A, at the Gardner transi¬ 
tion two-clones observables such as AAB{t) become dif¬ 
ferent from the long-time limit of dynamic observables 
such as A{t,tuj)- In this subsection, we thus estimate 
the location of the Gardner transition using an approach 
based on the study of two-clone static observables. These 
observables are static because, as we showed in Sec. IV A, 
they are time independent for p < pQ. We can thus 
arbitrarily choose any t to compute them. Here, we 
choose ts = 0.2V^/^ ^ 2, such that tq < G ^ T-metaj 
which we abbreviate below as A = A(G,tu, = 0) and 
Aab = AABits)- We use a similar notation for all other 


2. Probability distribution functions 

We next consider the probability distribution function 
(pdf) of the cage order parameters by computing A and 
Aab for each sample and constructing the histogram over 
samples. Note that because A is the mean square dis¬ 
placement restricted to a single sub-basin and the wait¬ 
ing time is fairly short, its distribution represents Psf(A) 
in the SF-RS phase, and only the peak around Aea in 
the SF-fullRSB phase (see Fig 2). 

Figure 10 shows P{A) and P{Aab) for po = 2.5 cal¬ 
culated from Ns = 40, 000 - 75,000 samples. The shape 
of P{A) is Gaussian-like at all p, and its mean value 
monotonically decreases with increasing p. The shape of 
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FIG. 9. (Color online). Density evolution of the averages 
of A and /S.ab for (po = 2.50. Numerical results (points) 
are compared with theoretical predictions for Ai, Aea and 
(A)sf (lines). Both datasets are scaled to reference values (see 
Sec. IV B 1 for details) to obtain a better agreement between 
theory and simulation. 


P{Aab), however, changes considerably over that same 
regime. For ip < pQ, it is Gaussian and analogous to 
that of f’(A), but near pQ it develops an exponential tail 
akin to a Gumbel distribution. If p is further increased, 
P{Aab) then becomes broader, which is consistent with 
the theoretical expectation (see Fig. 2). 

The development of an exponential tail at the critical 
point has been observed and studied for spin glasses in a 
field [59, 60]. The effect is thought to be due to disorder. 
Whereas the results for most samples fall within Gaus¬ 
sian fluctuations around a given mean value, a few rare 
samples have much larger Aab than the mean, giving 
an exponential tail to the distribution. The smaller the 
system, the stronger the effect (Fig. 11). These rare fluc¬ 
tuations are hypothesized to originate from the sample- 
to-sample fluctuations of the critical point [59], which 
then translates into significant sample-to-sample fluctu¬ 
ations of some of the measured observables. We come 
back to this point in Sec. IV G. 

The connection between the changing shape of the dis¬ 
tribution and criticality suggests that we can determine 
the critical transition from P{Aab) alone. We propose 
below two alternative procedures for detecting the Gard¬ 
ner transition using standard moments of the distribu¬ 
tion. 


3. Caging susceptibility 

We first define a caging susceptibility from the normal¬ 
ized variance of P[Aab) 


X = N 


(^ab) {^ab) 


(Aab) 


(24) 



0 0.01 0.02 0.03 0.04 0.05 0.06 

^AB 


FIG. 10. (Golor online), (a) P(A) and (b) P{Aab) at differ¬ 
ent p for po = 2.5. 


where the denominator corrects for the fact that (Aab) 
changes with p. As in the vicinity of any critical point, 
the susceptibility is expected to diverge as 

X(x{p^-p)-\ (25) 

where the critical exponent 1 is due to the fact that the 
MK model is mean-field in nature. 

The definition of x involves taking the quotient of two 
quantities that both suffer from strong finite-size correc¬ 
tions. In order to control for this effect we study the 
behavior of both terms as function of 1/N (Fig. 12). The 
denominator, (Aab), behaves as a regular series in 1/N, 
decreases smoothly with p, and eventually saturates 

above the Gardner point. The numerator, N{{A‘j^^) — 
-2 

(Aab) ) has, however, a more complex behavior. While 
it follows a nearly 1/N behavior for p < pQ with a 
small dependence on p, it grows significantly faster both 
with N and p for p > pQ. We attempt to extract the 
value of both quantities at the thermodynamic limit us¬ 
ing a second-order polynomial fit in 1/N, even though 

- - 2 

for p > pq the estimate oi N{{A/^^) — (Aab) ) obtained 
in this way is not reliable and larger systems would be 
needed to obtain a better extrapolation. 

The results for the susceptibility for N = 800 are re¬ 
ported in Fig. 13a for different values of po (and thus 
Pg)- For Po = 2.50 we also include x^^°° obtained 
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FIG. 11. (Color online), (a) Finite-size behavior of P{Aab) 
a.t ip = 3.0 « pa and (b) density evolution of the rescaled 
skewness F for three different N and for po = 2.5. 


using the extrapolations of Fig. 12. The comparison sug¬ 
gests that although the finite-size effects in both {Aab) 
- - 2 

and N{{A'\g) — (Aab) ) are still very strong for this sys¬ 
tem size, the determination of x is fairly well controlled, 
at least for p < pc- 

The numerical data in Fig. 13a suggest a roughly lin¬ 
ear behavior of x~^ for all po except po = 1.8, where 
the spacing between i^o and pc is narrowest, and where 
hopping most likely obfuscates the critical regime (see 
Fig. 16). For the other pq, we estimate the Gardner tran¬ 
sition by fitting Eq. (25) in the p < pc region (Table II 
and Fig. 16). 

In Fig. 14 we compare the numerical results with the 
theoretical predictions. Note, however, that the theoreti¬ 
cal curves do not correspond to the full x value, but only 
include the leading divergent term, namely the inverse 
of the replicon eigenvalue (see [16, 37] for details). As 
discussed for Fig. 9, both the theoretical and numerical 
datasets are also rescaled. The vertical axis is normal¬ 
ized to 1 at v? = (/JQ) while the horizontal p axis is scaled 
to Po, using the theoretical values in Table I and the fit 
values pj, for the numerics (Table II). Interestingly, the 
theoretical curves are not linear over the whole density 
regime. They instead bend towards po^ before vanish¬ 
ing linearly, but only in a fairly small region around pQ. 
The linear scaling observed in the simulation data is thus 
likely due to finite-size effects. The numerical estimates 
of Pq obtained from the linear fit of the whole set of data 
for p < PQ thus slightly overestimate the true transition 
point. Yet the effect is quite small - approximately 2%, 
based on Fig. 14. For p ^ pQ and above, we observe that 
finite-size effects dominate the numerical determination 


FIG. 12. (Golor online). Finite-size behavior of the two terms 
of the definition of the susceptibility in Eq. (24): (a) (Aab), 

and (b) N{{A\g) — (Aab) ), for po = 2.50. The dashed 
lines are fits of a second-order polynomial in 1/N to the data. 
The 1/N —>■ 0 extrapolation of these two fits have been used 
to extract the A —>■ oo limit of y, given in Fig. 13 by a dark 
solid line. 


of X, which remains finite instead of diverging. 


4- Caging skewness 

Near the Gardner transition, large sample-to-sample 
fluctuations give rise to a strong exponential tail in 
P{Aab)- This effect can be quantified by the skewness 
of the distribution 


^ (^^AB — (^ab)'J 

Recall that the skewness is a measure of a distribution’s 
asymmetry and that a Gaussian distribution would have 
F = 0. Sample-to-sample fluctuations are expected to 
be maximal at the critical point (see Sec. IV C), which 
provides an estimate of the Gardner transition, (see 
Fig. 13b and Table II). For all po, we see that is very 
close to the fitted divergence of the susceptibility, Pq. 
Finite-size analysis further shows that the rescaled skew¬ 
ness, r-\/]V, collapses the data for different N (Fig. 11). 
The peak of F is thus expected to persist all the way to 
the thermodynamic limit, consistently with comparable 
observations in spin-glass models [59]. 


F = 


(^Aab — (Aab)J 
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FIG. 13. (Color online), (a) Inverse caging susceptibility, 
X~^, and (b) skewness, F, as functions of ip for different po. 
Upon approaching pc, x apparently diverges as |</3 — 
except for po = 1.8. Note that the divergence of x coincides 
with the maximum of F. The N ^ oo extrapolation of the 
curves in Fig. 12 were used to obtain the dotted line. 

C. Sample-to-sample fluctuations 

We have assumed above that the abnormal behavior 
of r around the Gardner transition is due to sample- 
to-sample fluctuations, and we further motivate this hy¬ 
pothesis here. Recall that in our notation a sample is 
defined by a given configuration {ri} and a set of random 
shifts {A.,} (see Sec. Ill A). To compute the moments 
X and r of each sample, we perform Ng = 10, 000 inde¬ 
pendent clonings for each sample. In Fig. 15, the data 
for the ensemble of samples (same data as in Fig. 13) are 
compared with those of four individual samples. We note 
that the density evolution of F for the individual samples 
can have a very different behavior from the ensemble one. 
In particular, a peak around tpo is generally not seen. 

Sample-to-sample fluctuations have a smaller effect on 
X (Fig. 15). While the magnitude and the critical density 
exhibit some fluctuations, they all display a divergence 
similar to that of Eq. (25). Because each realization 



FIG. 14. (Color online). The numerical susceptibilities (scat¬ 
ters) are compared with theoretical predictions (lines). Both 
datasets are scaled to reference values (see Sec. IV B 3 for de¬ 
tails) to obtain a better agreement between theory and sim¬ 
ulation. We also report the linear fits to the numerical data 
used to determine Pq (dashed lines). This linear fit overesti¬ 
mates the transition point by roughly 2%. 



2.6 2.7 2.8 2.9 3 3.1 

<P 

FIG. 15. (Color online), (a) The susceptibility x and (b) the 
skewness F averaged over the ensemble of clones are compared 
with those of fonr individual samples, for ipo = 2.50. 


of disorder corresponds to different SF-RS metabasins, 
our results suggest that the metabasins themselves have 
slightly different properties. In particular, they exhibit 
different ipc, which is likely the physical origin of the 
exponential tail of P(Aab) and thus of its anomalous 
skewness [59]. 















15 



FIG. 16. (Color online). Inverse reduced pressure 1/p density 
ip phase diagram of the MK model in d = 3. The liquid EOS 
follows Eq. (6), while a specific state follows a glass EOS 
from ipo up to jamming at pj. Numerical estimates of the 
Gardner transition evolve with ipo similarly as the theoretical 
predictions [17]. 


V. SUMMARY OF RESULTS 

The Gardner transition at (pc was independently and 
quantitatively identified from: (i) the power-law diver¬ 
gence of the characteristic time r in Fig. 7, (ii) the linear 
vanishing of the inverse susceptibility x~^ iii Fig. 13a, 
and (Hi) the maximum of the skewness F in Fig. 13b. 
Table II contains these results and Figure 16 presents 
them in a phase diagram. The different estimates of pQ 
are generally quantitatively consistent with each other 
and qualitatively consistent with the mean-field SF re¬ 
sults from d ^ GO. 

The agreement with the mean-field calculation im¬ 
proves with density, likely because the Gaussian caging 
approximation used in the theory becomes a better ap¬ 
proximation at higher densities [34, 48]. A couple of 
reasons underlie the discrepancy between numerical es¬ 
timates and theory in the vicinity of pd- First, caging 
becomes imperfect at such low densities, which allows 
particles to hop between neighboring cages on a timescale 
comparable with the simulation time [34]. Hopping thus 
affects the dynamics of the system (A(t, tw) and AAB{t)) 
and also transforms the transitions (at either pd or pq) in 
crossovers. Second, because pc is expected to converge 
to pd upon approaching the dynamical glass transition 
(see Figs. 1 and 16), the critical regime becomes too small 
to make any fit to a critical power-law scaling. 

In addition to (i)-(iii), a more indirect way (iv) to de¬ 
termine Pq has been reported in Fig. 9 by comparing the 
density evolution of the cage order parameters with the 
theoretical predictions in d —>■ oo. However, while the 
effect is qualitatively (and visually) clear, obtaining a 
quantitative measure of pQ from this approach is some¬ 
what ill-defined. For p > pQ one also should treat A 


TABLE II. Results for pa and pj for various po- Errors 
are estimated as follows. For (/Jq, the error is the average 
distance between the susceptibility maximum and the next 
largest point. For Pq and Pq, fitting p intervals are varied to 
obtain the smallest and largest pa for which a fit is possible. 
The distance between the two values is the error. 


Pq 

r 

‘Pg 


‘Pg 

T 

PQ 

1.8 

2.10(5) 

2.11(3) 

- 

- 

2.534(3) 

2.2 

2.64(5) 

2.670(11) 

2.72(4) 

- 

2.876(4) 

2.5 

2.995(15) 

3.006(6) 

3.06(2) 

2.96 

3.151(3) 

3.0 

3.54(2) 

3.537(3) 

3.554(17) 

3.49 

3.622(5) 

4.0 

4.550(12) 

4.551(5) 

4.57(2) 

4.48 

4.584(4) 


results with caution. As discussed above, because of the 
finite-size and out-of-equilibrium nature of the system, 
the value of A drifts with time as different sub-basins 
are explored. Because in practice we evaluate A at a 
relatively short time 4 (Sec. IV B), we likely obtain a 
reasonable estimate of the size of a single sub-basin, but 
this procedure is also somewhat uncontrolled. 

Overall, we find that the most reliable way to deter¬ 
mine pG is the divergence of y (procedure (i)). This ef¬ 
fect is clearly the most spectacular signature of the tran¬ 
sition, and sample-to-sample fluctuations do not much 
affect its detection (Fig. 15). Interestingly, because xab 
is almost independent of time, one can determine y reli¬ 
ably using Aab at short times, as we did in this paper. 
Once PQ is determined in this way, a useful test is to 
check that this value is consistent with the behavior of 
A and Aab, as in Fig. 9. Procedure (Hi), i.e. finding 
the maximum of the skewness, requires averaging over 
a large number of samples which will surely be difficult 
in numerical simulations of realistic models of glasses as 
well as in experiments, where producing equilibrium con¬ 
figurations is extremely difficult. Procedure (i), i.e. the 
divergence of Tp, is also difficult because the determina¬ 
tion of T/d is subject to some ambiguity, but the study of 
the dynamics is useful because aging effects are manifest 
for p > PG (Fig. 5). Procedure (iv) is used here as a 
consistency check with the theory, rather than a method 
to detect pQ- 


VI. CONCLUSIONS 

The Gardner transition separates a stable glass 
metabasin at low densities (high temperatures) from a 
complex hierarchy of marginally stable sub-basins, at 
high densities (low temperatures). Its existence, which is 
proven in mean-field glass models [18, 19], has deep con¬ 
sequences on the low-temperature physics of glasses and 
on jamming [15, 26]. It is therefore extremely interest¬ 
ing to check whether such a transition exists in realistic 
models of glasses and in experiments. 

In this work we have investigated and compared several 
numerical procedures for detecting the Gardner transi- 
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tion in the MK model, which belongs to the universal¬ 
ity class of mean-field spin glasses while remaining fairly 
close to realistic glass formers [32]. We have presented 
three independent approaches for locating the transition, 
all of which show that the transition exists and is found in 
a region that is roughly consistent with theoretical pre¬ 
dictions. We have discussed the advantages and draw¬ 
backs of each of these strategies and the importance of 
finite-size effects. 

This work paves the way for studying the Gardner 
transition in more realistic numerical models of glasses, 
where the very existence of the Gardner transition is de¬ 
bated [29]. Our approach is also suitable to be repro¬ 
duced in experiments. SF, for instance, corresponds to 
a straightforward annealing, and some of the observables 
should be readily available through standard microscopy 
or scattering techniques. 

One key hurdle to generalizing our methodology to 
other systems is the need to equilibrate, without plant¬ 
ing (and thus through slow annealing), a glass state well 
above the (avoided) dynamic glass transition. To follow 
adiabatically a glass state, one should be able to prepare 
initial configurations such that the a-relaxation time 
is very large (tq, ^ rg), so that the SF experiment can 
be performed on time scale T/j ^ r ^ Tq,, as discussed 
in Sec. IID. For numerical simulations this requirement 
can be particularly computationally onerous, but it may 
be more easily achievable in experimental systems, where 
longer timescales can typically be reached. In particular, 
it would be very interesting to investigate the existence 
of the Gardner transition in ultrastable glasses that can 
be prepared through vapor deposition and have an ex¬ 
tremely large Tq, [51, 52, 61, 62]. In experiments, the 
bigger challenge would be to substitute the cloning pro¬ 


cedure with a (potentially very) large number of experi¬ 
mental replicates. 

Finite-dimensional non-mean-field glass formers dis¬ 
play features that are not observable in the MK model. 
In particular, we expect a diverging length scale to be 
associated with the Gardner transition in these systems. 
This length scale is expected to capture static hetero¬ 
geneity, which represents the spatial inhomogeneity of 
cage sizes around and above (^g- In principle, this kind 
of static heterogeneity should be different from both the 
dynamic heterogeneity around the dynamic glass transi¬ 
tion, and the heterogeneity close to jamming, which is 
related to soft relaxation modes [63]. Understanding the 
relevance of marginal stability for glassy dynamics, and 
the relation with the ideas of Ref. [63], would open a new 
window on the property of low-temperature glasses. 
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